Josephson current through a single Anderson impurity coupled to BCS leads 
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We investigate the Josephson current (J(<?!>)) tlirougli a quantum dot embedded between two 
superconductors showing a phase difference 0. The system is modeled as a single Anderson impurity 
coupled to BCS leads, and the functional and the numerical renormalization group frameworks are 
employed to treat the local Coulomb interaction U. We reestablish the picture of a quantum phase 
transition occurring if the ratio between the Kondo temperature Tk and the superconducting energy 
gap A or, at appropriate Tk/A, the phase difference (j) or the impurity energy is varied. We present 
accurate zero- as well as finite-temperature T data for the current itself, thereby settling a dispute 
raised about its magnitude. For small to intermediate U and at T = the truncated functional 
renormalization group is demonstrated to produce reliable results without the need to implement 
demanding numerics. It thus provides a tool to extract characteristics from experimental current- 
voltage measurements. 

PACS numbers: 74.50.-|-r, 75.20.Hr 
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I. INTRODUCTION 

The Kondo efFcct^'^ and superconductivity"^ arc two of 
the most striking manifestations of electronic correlations 
in low-energy condensed matter physics. The interplay 
of both phenomena, e.g. showing up for superconducting 
metals containing magnetic impurities, was first studied 
several decades ago.^"^^ The Kondo temperature Tk and 
the superconducting gap A are the two competing energy 
scales which govern the low-energy physics of such sys- 
tems. If Tji ^ A, local magnetic moments are screened 
by virtue of the Kondo effect. This causes Cooper pairs 
to break, and the ground state becomes a Kondo rather 
than a BCS singlet. In the opposite limit Tk ^ A, 
Kondo screening is disturbed due to the superconduct- 
ing gap at the Fermi energy. The ground state describes 
free magnetic moments. At temperature T = 0, a quan- 
tum phase transition from a non-magnetic singlet to a 
degenerate (so-called magnetic) ground state is observed 
if A/Tk increases. 

In recent years, a renewed theoretical interest in the in- 
terplay between Kondo and BCS physics has developed 
due to the rise of nanotcchnology and the associated re- 
alization of quantum dot systems connected to supercon- 
ducting leads. ^^^^"^ The microscopic parameters of such 
nanoscale systems (e.g. the energy e of the quantum dot) 
can be easily tuned, thereby allowing to study the physics 
in a controlled way. In general, an equilibrium current 
( J((/))) flows through such a quantum dot Josephson junc- 
tion, provided that there is a finite phase difference (fi be- 
tween both superconductors. From the theoretical point 
of view, the single impurity Anderson model with BCS 
superconducting leads can be used to describe the low- 
energy physics. If the local interaction U between spin 
up and down electrons is sufficiently large so that the 
impurity is singly occupied, there is again a competition 
between Kondo screening and the formation of Cooper 
pairs. The Kondo singlet ground state becomes a mag- 



netic doublet if A/Tk is increased at arbitrary phase 
difference (f>. However, the critical value Uc{A) describ- 
ing the phase boundary depends on cf). Hence, a phase 
transition can be observed if the phase difference is var- 
ied at appropriate A/Tk- Likewise, a transition to the 
singlet state is observed if the system is driven away 
from particle-hole symmetry by a gate voltage e. At the 
critical point, the sign of (J) changes discontinuously at 
T = 0. 

If one takes the limit A — > 0, the single impurity 
Anderson model with ordinary Fermi liquid leads is re- 
covered. The latter is well-known to describe strongly- 
correlated electron behavior. Hence, reliable many- 
particle methods are needed for an accurate treatment 
of the interaction U between electrons occupying the 
impurity. In this paper, the functional (fRG) and nu- 
merical (NRG) renormalization group schemes will be 
employed. We mainly focus on the zero-temperature 
limit but towards the end of the paper also study fi- 
nite temperature effects. The NRG is a very reliable 
method to investigate physical properties of systems with 
Coulomb interaction^^ and thus provides a powerful tool 
for an unbiased calculation of the Josephson current. 
Unfortunately, it requires large numerical effort and in 
practice only small systems of high symmetry can be 
treated. In contrast, truncated fRG is fast, flexible, easy- 
to-implemcnt, and free of numerical parameters, but by 
construction limited to small to intermediate interaction 
strength. Comparison with NRG data, however, showed 
that fRG correctly describes zero-temperature (i.e. zero- 
frequency) transport properties of multi-level quantum 
dot geometries connected to Fermi liquid leads up to 
fairly large [/.^^'^^ In this paper we establish the accuracy 
of fRG in treating the Josephson problem by comparing 
with analytical results at A = oo as well as with our 
NRG data. 

Experimentally, it is impossible to control the phase 
difference (p so that only the weighted one-period aver- 
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age of the Josephson current can be measured. In or- 
der to extract physical properties, the experimental data 
needs to be fitted by a procedure which requires current- 
phase characteristics as an input parameter. This paper 
shows that the fRG is an accurate and fast tool to provide 
(J(0)). 

As mentioned above, the fundamental physics of mag- 
netic impurities inside a superconductor was explained 
decades ago.'*"^'^ Due to the revived interest motivated 
by recent experiments, ^^■^^'^^•^'^ the quantum dot Joseph- 
son junction has been studied using various theoretical 
approaches. In particular, the NRG was applied to ac- 
curately calculate the phase boundary between the sin- 
glet and doublet phase^*'^^ as well as the single-particle 
spectral bmction.'^'' Even though the phase transition 
is already captured on the Hartrcc-Fock IcveP"'^'^^ and 
by pcrturbative approaches, "^^'^^ there are few quanti- 
tatively reliable results for the Josephson current at 
T = 0. The atomic limit A = oo can be treated ana- 
lytically (see Sec. VB), and Glazman and Matveev de- 
rived an expression for ( J(</>)) in the limits of F — > and 
A 0, respectively.^^ Quantum Monte Carlo (QMC) 
calculations were carried out by Siano and Egger (SE, 
Ref. 36), but they show significant finite-temperature 
effects (QMC being an inherently finite-T - method). 
Sellier et al. published finitc-tcmpcrature data for the 
infinite-?/ Anderson model obtained from the noncrossing 
approximation."^^ NRG results for {J{(t>)) at arbitrary pa- 
rameters were published by Choi et al. (CLKB, Ref. 38), 
but they have been criticized by SE to be inaccurate. 
The present paper settles this dispute. 

We organize this paper as follows. In Sec. II we intro- 
duce our model Hamiltonian. In Sec. Ill we briefly com- 
ment on the general ideas of the fRG and NRG. In par- 
ticular, we derive the fRG flow equations for the present 
context. We present our results for the Josephson current 
in Sec. IV. Sec. V is devoted to an extended discussion 
of different theoretical approaches to the Josephson prob- 
lem. In particular, we check our NRG numerics against 
the exact solution at J7 = 0. Analytic results for A = oo 
provide a benchmark for the fRG approximation. We 
comment on the dispute between CLKB and SE. In the 
appendices we present details of the NRG calculation, de- 
rive an exact formula for the Josephson current in terms 
of the self-energy, and comment on the issue of current 
conservation. 



II. THE MODEL 

As indicated in Fig. 1 , our model Hamiltonian consists 
of three parts: 



coup 



(1) 



s=L,R 



s=L.R 



The first term describes a single Anderson impurity (the 
quantum dot) with onsite energy e and Coulomb rcpul- 
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FIG. 1: (Color online) The quantum dot Josephson junction 
considered in this paper. 



sion U between spin up and spin down electrons: 



dot rjdot 



= h: 



dot 



(2) 



The single-particle energy was shifted such that e = 
corresponds to the point of particle-hole symmetry. The 
left (s = L) and right (s = R) superconducting leads are 
modeled by a standard BCS Hamiltonian, 



#s t t 



k(7 



H.c.j , 

(3) 

where As and (ps denote the magnitude and phase of the 
superconducting order parameter, and Cska are the 
usual annihilation operators of the dot and lead electrons, 
respectively. The quantum dot is coupled to the leads by 



(4) 



The hopping matrix clement tg is assumed to be real. We 
have introduced the local electron operator at the end of 
the leads, Csa = Y.k^ska / Vn . 

In order to derive the fRG flow equations, introduction 
of the Nambu formalism will prove to be helpful. To 
this end, we recast the Hamiltonian in terms of anti- 
commuting spinors. 



*,sfe = 



sfcT 



ay 
d\, 



(5) 



The dot part can then be written as 



H^°' = eV'Vg^ - U (^^1^1 - - , (6) 

with Lpi^2 denoting the components of the Nambu oper- 
ator If. We have introduced the Pauli matrices r^. For 
the BCS leads we obtain the usual expression 



with 



e*'^' 



(7) 



(8) 



-^4"^ 

Finally, the coupling Hamiltonian becomes 

^coup ^ ^tg^lr^^ + H.c. , (9) 
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where = sk/^/N . One should note that by in- 

troduction of the Nambu formaUsm, we get rid of all 
anomalous expectation values but have to deal with an 
additional single-particle quantum number (the Nambu 
index) . 



III. THE METHODS 

A. Functional RG 

The fRG is one implementation of the general renor- 
malization group idea for interacting quantum many- 
particle systems. It starts with introducing an energy 
cutoff A into the noninteracting Green function 5°. Here, 
we choose an infrared cutoff in Matsubara frequency 
space, 

G"'^{tuj) = e{\Lu\ - A)g'^{tuj). (10) 

Using this propagator, the m - particle vertex functions 
7^ acquire a A - dependence. Differentiating each 7^ 
with respect to A yields an infinite hierarchy of flow equa- 
tions 

5a7^. = fm iH}) . (11) 

The functions /,„ can be computed straight-forwardly by 
introducing a generating functional for 7,„. At A = 0, the 
original cutoff-free problem is recovered. Hence, one can 
obtain an exact expression for the sclf-encrgy S = —^i^^ 
by solving the set of coupled differential equations (11). 
In practice, this infinite hierarchy has to be truncated, 
thereby rendering the fRG an approximate method. In 
this paper, we will employ a truncation scheme that keeps 
the flow equations for the sclf-encrgy and the two-particle 
vertex 7^ evaluated at zero external frequency. The 
resulting approximation to the self-energy is frequency- 
independent and can be computed with minor numerical 
effort. By construction, this truncated fRG becomes ex- 
act in the noninteracting limit. One can show that it 
keeps track of all terms of first order in U but the RG 
procedure leads to an efficient resummation of higher or- 
der terms. The truncated fRG was demonstrated to suc- 
cessfully describe strong correlation effects (such as im- 
portant aspects of Kondo physics). 

At zero temperature, the approximate flow equations 
for := —ji and := 7^ explicitly read (a detailed 
derivation can be found in Ref. 42) 

ui = ±A 2,2' 



for the self-energy, and 

uj=±A 3,3' 4,4' I 
X r^',2';3.4r^',4';l,2 + 2^3^3, (zc^)^4^4, (zw) 

X ^ 2',4';l,3-'^ 3',1':4,2 ^ l',4';l,3^ 3',2';4,2 f 

(13) 

for the effective interaction. The lower indices denote 
arbitrary single-particle quantum numbers. We have de- 
fined 

[g^iiLo)]^' ^[g'^itco)]'' (14) 

One should note that the sharp cutoff has completely dis- 
appeared from the zero-temperature flow equations (12) 
and (13), rendering them easy to integrate numerically. 
In practice, it is convenient to exploit symmetries of the 
two-particle vertex (such as anti-symmetry and spin con- 
servation) in order to reduce the number of independent 
equations. 

Due to the static (frequency-independent) approxima- 
tion involved in setting up our fRG scheme one cannot 
expect to obtain reliable results at flnite temperatures 
for the problem at hand. Thus, we will focus exclu- 
sively on performing fRG calculations at T = 0, even 
though the flow equations (12) and (13) can be gener- 
alized straight-forwardly to the T > - case. Study- 
ing flnite-tcmperature effects with this method requires 
a more sophisticated truncation scheme and is subject of 
ongoing research. 

A numerical solution of the flow equations starts at 
some large but finite Aq. Due the slow decay of the r.h.s. 
of Eq. (12), the integration from A = 00 down to A = 
Aq always yields a finite contribution. This contribution 
tends to a constant if Aq — > 00 and can be computed 
analytically, leading to the following initial conditions at 
some large but finite Aq: 

SAn — *00 1 \ ^ - T^An — >00 — /'1 

1,1' = 2 Z^"l'.2;l,2 , ri,';2,.i^2 = «1',2';1,2, (15) 

2 

where v is the bare antisymmetrized two-particle inter- 
action. The convergence factor in Eq. (12) can then be 
dropped and the flow equations integrated numerically. 

Application of the fRG scheme to the quantum dot 
Josephson junction is achieved straight-forwardly within 
the Nambu formalism. The first step consists of integrat- 
ing out the noninteracting leads using a standard projec- 
tion technique. ^'^ Thereafter, instead of dealing with an 
infinite system we only have to consider two interact- 
ing (Nambu) particles. For the noninteracting dot Green 
function we obtain 

[g'^it^)] -'=ioj- ht°' - E isr35^(*^)^3, (16) 

s=L,_R 
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where /iq°* is the single-particle version of i?Q°* and g''{z) 
denotes the local propagator at the last site of the isolated 
BCS lead: 



lUJ 

-A,e- 



ILO 



(17) 



We assume the local density of states at the end of the su- 
perconducting leads Ps to be energy-independent (wide- 
band limit). The dot propagator then explicitly reads 



lOJ) 



where we have introduced 



A(ia;)* iu + e 



1 



E 



v4J^+Af 



(18) 



(19) 



with Ts = Trps\ts\'^ being the dot-lead hybridization, and 

r.A, 



A{iLu) 



E 



(20) 



In this paper we will mainly focus on the case of symmet- 
ric couplings (r^ = Tfj = r/2) and equal superconduct- 
ing gaps (Ai = Afl = A). If (without loss of generality) 
we additionally choose = — = 0/2, the function 
A{iuj) = rAcos((/)/2)/vuJ^-l-A^ becomes purely real. 
In general, the self-energy of the quantum dot Joseph- 
son junction is characterized by a diagonal component 
S S R and an off-diagonal (anomalous) part Ea G C. 
Within our truncated fRG approximation, they are both 
frequency- independent. Thus, the matrix reads 



1 I ILO - 



- A(icj) 



where we have defined the determinant 



D^{iio) ^Cj'^ + {e + E^) + A{iuj) 



iCo - 



A{iu}) 
e-E^ 



(21) 



(22) 



The flow equations (12) for the self-energy can now be 
written explicitly as 



SaE^ 



for the diagonal component, and 

Ei - A{ih) 



(23) 



(24) 



for the anomalous part. In the symmetric case, E^ is 
real (since A is real). denotes the only independent 
component of the two-particle vertex. Its flow equation 
(13) is given by 





2 p 







A{iK) 



(25) 



An even simpler approximation can be obtained by ne- 
glecting the flow of the two-particle vertex altogether, 
using a constant interaction ~ U . While this does 
not affect our results qualitatively, the accuracy of fRG 
significantly improves if the flow of is accounted for 
(see Sec. IV B). The initial conditions read 



■f^Ao^oo 
^A 



, u 



Aq— >c>o 



u 



(26) 



In order to obtain the fRG approximation to the self- 
energy E = E'^^^ and Ea = ^a^^j we solve these 
coupled differential equations by numerically integrating 
from Aq = 10^ down to A = using standard routines. 

Having obtained the self-energy of the system, the 
Josephson current at T = can be computed from the 
following integral (here focusing on A^, = A^j = A, 
Tl = r_R, and 0^ = -0_r = 0/2): 



r2A2sin(0) 2rAEAsin(0/2) 



D{iLu) 



du>, 
(27) 



where D{iuj) — D^^*^{iuj). We take h — 1 and e — 1 
(the latter being the electron charge) in the following. 
We will derive Eq. (27) and its generalization for T > 0, 
non-symmetric leads and a self-energy which is explicitly 
frequency-dependent in Appendix B. In general, cur- 
rent conservation (J^) — —{Jr) is ensured for A^ = Ar 
(and otherwise arbitrary parameters) within our fRG ap- 
proximation. The general issue of current conservation is 
commented on in Appendix C. 



B. Numerical RG 

The main idea of the NRG in application to quan- 
tum impurity systems is to discretize the flat conduc- 
tion band of the leads using a set of logarithmic energies 
{±£)A~", n > 0}, with D being half the bandwidth and 
A > 1 the NRG discretization parameter. ^•^"'■^^ The re- 
sulting discrete model can be mapped onto a semi-infinite 
tight-binding chain with the impurity being the first site. 
The Hamiltonian of the semi-infinite chain is then diago- 
nalized iteratively by adding one site at a time, starting 
out with the isolated impurity. Due to the logarithmic 
discretization, the hopping matrix elements u„ between 
successive sites fall off exponentially with the distance n 
from the impurity (u„ ~ A~"/^), rendering it possible 
to resolve smaller and smaller energy scales during the 
iteration. For a numerical implementation, a truncation 
procedure has to be employed as the dimension of the 
Hilbert space grows exponentially with the length of the 
chain. A very simple truncation scheme is to retain only 
the A^c lowest-lying many-particle states at each iterative 
step, which is reasonable since the states of the shorter 
chain affect those of the longer only in a small energy 
window ~ A^^/^. 

The essential approximations in the NRG approach 
arc the logarithmic discretization of the conduction band 
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and the truncation of the Hilbert space during the it- 
erative diagonalization, implying that the accuracy of 
this method is controlled by the parameters A and Nc- 
For models containing only a single conduction band (so- 
called single-channel models) , NRG calculations at A = 2 
and Nc = 500 typically agree well with analytical (such 
as Bethe ansatz) results. In this paper, we are consider- 
ing a two-channel model which could cause problems as 
the number of low-energy states increases exponentially 
with the number of channels. Hence, NRG calculations 
have to be performed with care and checks at different 
A and Nc against analytical results (available at [/ = 0) 
are imperative (see Sec. VA). 

The NRG in the present work was implemented such 
that two sites from different channels were simultane- 
ously added at each iteration, and we have kept the lowest 
Nc = 700 energy states after the diagonalization proce- 
dure. Then the dimension of the Hamiltonian matrix in 
each recursive NRG step becomes approximately 700x4^. 
The effective dimension in actual calculations has been 
reduced efficiently by the use of symmetries, which is 
helpful for improving the numerical accuracy. Namely, 
in the case of equal gaps (A^ = A^) and symmetric cou- 
plings (Fi = Ffl), the Hamiltonian defined in Eq. (1) can 
be described by a real symmetric matrix even if the phase 
difference (j) is finite (see Appendix A). Furthermore, in 
the particle- hole symmetric case (e = 0), there is an addi- 
tional U(l) symmetry associated with the conservation of 
the pseudo-spin variable Q^^ [Eq. (A8)]. We have carried 
out our NRG calculations using the quantum numbers Q 
and 5". The former is related to Q^^ whereas the latter 
corresponds to the SU(2) symmetry of the real spin. We 
have checked whether Nc = 700 is enough for obtaining 
accurate results by increasing the number of kept states 
up to Nc = 1716. It turned out that in contrast to single- 
channel models in our case it is not sufficient to retain 
only Nc = 700 states at A = 2, particularly for F > A. 
We have also confirmed that ground-state properties can 
be obtained reasonably well with Nc = 700 kept states 
for A > 4 (for the details see Sec. V A). The addition of a 
new site during the iteration procedure can be viewed as 
a perturbation of relative strength of the order of A~^/^ 
(specifically for normal leads A = 0) , implying that with 
decreasing A one has to keep more and more states to 
get reliable results. Hence, the numerical accuracy can 
be improved by increasing A at fixed Nc, even though 
the approximation involved becomes exact in the oppo- 
site limit A ^ 1. For F < A, the wave function of the 
Andreev bound state is localized well at the impurity 
site and does not penetrate deep inside the BCS leads. 
In such cases the convergence becomes better, and Nc 
can be smaller than the value one needs for achieving the 
same accuracy in the opposite case F ^ A. We note 
that in our calculations the NRG Hamiltonian given in 
Eq. (A2) is diagonalized exactly (without introducing the 
truncation) up to 7 sites, which consist of the dot at the 
center, and the first three orbitals (n = 0, 1, 2) of the s„ 
and Pn particles, respectively (see Appendix A). There- 
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FIG. 2: (Color online) The phase diagram at e = 0, = 0, 
and Fl = characterized by tlie boundary line between 
the singlet (upper region) and doublet (lower region) phase. 
Solid (dashed) lines show fRG results with (without) flow of 
the static vertex obtained from Eqs. (23), (24), (25), and (27) 
and the definition of the phases via J(<i> — > 0) ^ 0. Sym- 
bols are NRG data (taken from Refs. 28 and 29). The inset 
shows that near the origin (at large A), the slope of the phase 
boundary is 1/2 (dashed line) in accordance with analytical 
results obtained in the atomic limit [see Eq. (38)]. 



fore, the truncation only starts at the next NRG step for 
the cluster with 9 sites. 

For the NRG calculations shown in Fig. 4(a) and 
Fig. 6, we have chosen the half bandwidth such that 
F = 0.03847D. We have used a different parameter for 
Fig. 4(b), namely F = 0.0021?. The flow of low-energy 
eigenvalues converges after at most 50 NRG steps, and 
thus the iterations up to n < 50 were enough for the 
parameter sets we have examined. The current was com- 
puted as the expectation value of the discretized current 
operator given in Eq. (A7). 



IV. NUMERICAL RESULTS AT T = 

In this section, we present our numerical results. 
Within fRG, the singlet and doublet phase is defined 
via (J) > and (J) < (focusing on < < tt), 
respectively. Within the NRG framework, the degener- 
acy of the ground state is directly accessible. First, we 
report on what happens during the integration of the 
fRG flow equations. In particular, we show that the cur- 
rent is fairly sinusoidal in the doublet phase. Next, we 
show phase diagrams at (e = 0) and away from (e ^ 0) 
particle-hole symmetry which we obtain using fRG and 
NRG. The NRG data was previously published by one of 
us (A.O., Refs. 28 and 29). We stick to equal supercon- 
ducting gaps (Al = Afl = A) but study both symmetric 
(Fi = Ffl) and asymmetric (F^ 7^ Vb) couplings. As 
all physical quantities (such as (J)) depend on the phase 
difference only, we choose 0^ = —4'R = 0/2 from now on. 
We present new fRG and NRG results for the Josephson 
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FIG. 3: (Color online) The same as in Fig. 2, but for both symmetric Tl ~ Tr (dashed lines) and asymmetric Fl ~ 1.44rH (solid 
lines) at (a) ^ and (b) — > tt. Symbols denote NRG data. The dotted lines are the phase boundaries for e — — 2.5A + U/2, 
Tl = IAATr, and <?!> = tt. 



current J(0) = ( J((/>)) and demonstrate that for weak to 
intermediate interactions fRG performs well in compar- 
ison with NRG in describing both the phase boundary 
and the magnitude of J. 



A. Integrating the fRG flow equations 

First, we give a brief overview of what is happening 
during the integration of the fRG flow equations. For 
simplicity, we stick to F^ = F/j = r/2. In the singlet 
phase, nothing particular happens when (23), (24), and 
(25) are solved numerically. In the doublet phase, one 
observes that at a certain Ac, e + T.^" becomes zero, 
implying that T,^ = -e for all A < Ac [see Eq. (23)]. 
The anomalous component continues to flow and reaches 
Ea = Fcos((/)/2) at A = 0. Hence, the doublet phase can 
already be identified during the flow. Plugging S = — e 
and Ea = Fcos(^/2) into Eq. (27) yields an analytic ex- 
pression for the Josephson current in the doublet phase: 



J = 



F2Asin(0) 



A 



1 



A2 



duj, (28) 



which is a universal curve independent of both U and e. 
One should note that sin(0) is not the only (jj - depen- 
dence as 



Ddiito) ^lD^ cos2((/)/2) 1 



(29) 



Even though the complete independence of U and e is 
an artifact of the fRG approximation, it is instructive 
to gain analytical insight into the current described by 
Eq. (28). To this end, we scale A out of the integrand: 



J 



-FV(27rA)sin(0) {y - y^) 



x2(l + y(F/A))2 + (F/A)2 cos2(</)/2)(l - y)- 



2 dx, 
(30) 



with y = 1/Vl + x"^ 
behavior at small \x 



This integral is dominated by the 
The term proportional to cos((/)/2) 
in the denominator is then of order and can be ne- 
glected compared to the term, provided that F/A is 
not too large (which is generally the case in the doublet 
phase). This gives 



lim J 

r/A^o 



27rA 



1 



1 



A 271 



■ sin((/)) 



x'^ 1 + ; 
F2 

-0.18— sin(0). 



dx 



(31) 



Thus, the current in the doublet phase obtained from fRG 
is fairly sinusoidal (becoming perfectly sinusoidal with 
decreasing F/A) and of order F^/A (this is illustrated 
by Fig. 8; see below). 

Within fRG, we compute the average occupation num- 
ber Ud of the quantum dot from integrating the Green 
function Qi^i — Gif^ over the imaginary axis: 



nd 



^ I e"^'>gi,iiiu^)dLU. 
Ztt J 



(32) 



In the doublet phase, t/i,i(ia;) becomes an odd function 
(since E + e = 0), implying that n = 1/2 due to the 
contribution from large frequencies. It is again an arti- 
fact of the truncated fRG that in this phase the average 
occupation is completely independent of U and e. How- 
ever, NRG computations showed that the deviations from 
ridie) = 1/2 arc small. 



B. Phase diagrams 

Fig. 2 shows the phase diagram at the point of particle- 
hole symmetry (e = 0) for zero phase difference and sym- 
metric couplings. It confirms that if cither U is increased 
or F decreased, implying that the Bethe ansatz Kondo 
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FIG. 4: (Color online) Zero-temperature Josephson current J in units of Jc = eA/ft as a function of the phase difference 4> 
computed with fRG (solid lines) and NRG (symbols) at e = and Tz, = r^. (a) A is varied at fixed U /V = 5.2 [Tk/V = 0.209). 
The parameters correspond to = 0.11, = 1.76, and A/Tk ~ 11 and were chosen for direct comparison with Fig. 3 

of Ref. 38 (please note Ref. 46). For clarity, the curves at A/Tk = 11 were scaled up by a factor of 20. (b) A/F — 0.5 is fixed 
at different A/Tk ~ 1.1, A/Tk = 1.7, and A/Tk = 5.8. The discretization parameter for NRG was chosen to be A = 4 (for 
details see Sees. TUB, VA, and Appendix A). 



scale 



Tk - Vc7r72cxp 



8ur 



I [/2 -46^1 



(33) 



decreases, the system shows a phase transition from 
a non-magnetic singlet to a magnetic doublet ground 
state.*'' If U and T are fixed, the system is always in a sin- 
glet state at sufficiently small A. As A becomes larger, 
a phase transition is observed provided that U > 2r. 
In contrast, for T > U/2 the Kondo effect is not ac- 
tive and the ground state always remains a (BCS) singlet 
no matter how large A. This can be understood in de- 
tail from the analytical treatment of the so-called atomic 
limit A = oo. The inset of Fig. 2 illustrates that at large 
A, the phase boundary indeed approaches the analytical 
result r^iU) = U/2 [see Eq. (38) of Sec. VB]. 

A phase diagram for and finite phase dif- 

ference (j) is shown in Fig. 3. The general picture of the 
phase transition remains the same, only the phase bound- 
ary is affected. One can see that asymmetric couplings 
Tl 7^ Ffl stabilize the singlet phase. In contrast, the 
doublet phase becomes more and more favorable with 
increasing (j). The phase boundary continuously evolves 
from the = into the cf) = n curve if (f) is gradually 
increased from to tt. If F^ = F^ and (/> = tt, the sin- 
glet phase completely disappears (both within £RG and 
NRG). The phase boundary at finite e (dotted lines) ap- 
proximately starts from U = 2.5 A, which is the point 
where in the atomic limit (F/A ~ 0) singlet and dou- 
blet states cross. Generalizing the arguments presented 
in Sec. VB for F^ 7^ ^r,^^ it is possible to demonstrate 
that the boundary line shows square-root behaviour close 
to F = 0. One should note that in our case the deviation 
of the dot energy from the point of particle-hole sym- 
metry depends on the interaction strength. Hence, the 
Kondo temperature Tx given by Eq. (33) increases with 



U , causing the non-monotonicity in the phase boundary. 
Additional NRG results away from particle-hole symme- 
try can be found in Ref. 29. 

In general, the phase boundaries obtained from fRG 
and NRG agree quite well. As expected, the agreement 
is particularly good at small U. The dashed line in the 
main part of Fig. 2 shows fRG results where the flow 
of the static vertex was discarded. Whereas the general 
physics is captured by this even simpler fRG truncation 
scheme (which sets = U), the quantitative agreement 
with NRG improves if the flow equation (25) is accounted 
for. Hence, we will stick to this improved approximation 
from now on. 



C. Josephson current 

In Fig. 4, we show fRG and NRG results for the Joseph- 
son current J as a function of the phase difference (j). The 
current is given in units of Jc — eA/h. The flgure illus- 
trates the generic physical picture; in particular, we ob- 
serve the same away from particle- hole symmetry (e 7^ 0) 
and for asymmetric couplings (F^ 7^ Fjj). In certain lim- 
its analytical results help to understand the form of J{4>). 
We will discuss this in detail in Sec. V. However, it is 
instructive to recall that the current flowing between two 
superconductors connected by a weak tunneling barrier 
is purely sinusoidal [J ^ sin((/))].*^ 

As discussed above, the doublet phase is stabilized if 
U, A, or (j) is increased. Thus, at appropriate A/Tk, the 
phase transition manifests as a discontinuity at 4>c in the 
J{<j)) curve. This is illustrated in Fig. 4(a,b) (at fixed 
Tk and A, respectively). Note that the parameters of 
Fig. 4(a) were chosen such that we can directly compare 
our results with those of Ref. 38 (see Sec. V C).*^ As there 
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is a difference between the phase boundaries obtained 
from fRG and NRG, (l>c is different in both approaches. 
However, the form of the curves is captured well by fRG 
even at large U/T. The agreement between fRG and 
NRG concerning the current amplitude is good at small 
to intermediate U, being almost perfect throughout the 
singlet phase. This is satisfying as our truncated fRG is 
by construction well-controlled in this regime. The agree- 
ment becomes worse as U becomes large [see Fig. 4(b)]. 
In general, the deviation between both methods is most 
severe in the doublet phase at intermediate A, becoming 
better as A increases [compare Figs. 4(a) and (b)]. It 
is an artifact of the fRG approximation that the current 
is completely independent of U. However, the sinusoidal 
form of J{<f>) in the doublet phase is described equally 
well by both methods. 

One parameter that can be easily tuned experimen- 
tally is the energy of the quantum dot,^^ rendering it 
reasonable to consider the current as a function of e. 
Within fRG, computing J{4>) away from half- filling is 
not different from computing J{4>) at e = 0. In con- 
trast, NRG greatly benefits from symmetry properties 
which only hold at e = (see Appendix A) , and comput- 
ing the current away from this point is (though possible) 
numerically demanding. As we have demonstrated fRG 
to provide acceptable accuracy in comparison with NRG 
at small to intermediate f/, we refrain from this task and 
show only fRG results for J(e) and the average occupa- 
tion number n{e) in Fig. 5. As mentioned before, the 
singlet phase is always favored at large |e|. If A/T^"" 
is appropriately chosen, the phase transition to the dou- 
blet ground state occurs at a certain [ed- The system is 
in the doublet phase for |e| < [ed, and both the current 
and the occupation nd{e) = 1/2 are independent of U 
and e within the fRG approximation. NRG calculations 
showed that the deviations from nd{e.) = 1/2 are indeed 
small. One should note that our findings for J(e) are 
consistent with recent experiments.^^ 

Within fRG, another possibility to compute the 
Josephson current is to differentiate the grand canonical 
potential fl w.r.t. the phase difference cf) [see Eq. (C3)]. 
Since the truncated fRG is a non-conserving approxi- 
mation, the current computed in this way will generally 
differ from the one obtained via the self-energy formula 
Eq. (27). The quantity can be obtained directly from 
a flow equation. From the formalism one would expect 
the energy to be computed more accurately than the self- 
energy (the former being a vertex function of lower or- 
der), rendering it reasonable to calculate the Josephson 
current from the approximated energy rather than from 
the approximated self-energy. This expectation is contra- 
dicted by the observation that the former way does not 
give meaningful results in the doublet phase; in partic- 
ular, the current remains positive. In the singlet phase 
the current computed from the energy compares badly 
to exact results obtained at A — s- cx). Further studies on 
this issue are required. For the time being, we calculate 
J from the self-energy via Eq. (27). 
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FIG. 5: (Color online) fRG results for the Josephson current J 
and the average occupation number n at T = as a function 
of the dot's energy e at U/T = 3 (solid), U/T = 5 (long 
dashed) , and U/T = 7 (short dashed) . The other parameters 
are A/F = 1, (j)/n = 0.5, and Tl = Tr. These results are 
very much consistent with recent experiments (see Fig. 4 of 
Ref. 23). 



V. OTHER APPROACHES 

In this section we discuss different approaches to the 
Josephson problem. An analytical treatment is possible if 
either U = or A = oo. Both cases were partly adressed 
in prior works (see e.g. Refs. 29 and 33 for the atomic 
limit) but are again reported on here since the exactly- 
solvable noninteracting case is an important check for 
the NRG numerics, whereas analytical results at C/ > 
provide an additional benchmark for the fRG. We argue 
that a previous NRG approach by Choi et al. (Ref. 38) 
does not agree quantitatively with our results that we 
believe to be correct, although qualitatively both NRG 
approaches yield similar data except for a few obvious 
errors. We present NRG calculations of the Josephson 
current at finite temperatures and compare these results 
with the quantum Monte Carlo approach by Siano and 
Egger (Ref. 36). The discrepancy between both methods 
can be explained by considering the first excited many- 
body energies obtained from NRG. Finally, we comment 
on the issue of mean-field calculations. 



A. The noninteracting case 

At U = 0, the ground state of the system is always a 
(BCS) singlet. There is no doublet phase. An analytical 
expression for the zero-temperature Josephson current is 
provided by Eq. (27). Setting the self-energy to zero (and 
focusing on = Ffl), the integral can be rewritten as 

j = ± / ^^'"('^) dr 

2tt J cos2(0/2) + a-2(l + (A/r)y)V(e/r)2y2 ' 

(34) 
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FIG. 6: (Color online) NRG results for the zero-temperature 
Josephson current through a noninteracting, particle-hole 
symmetric dot for A/F = 0.023076 and several sets of the 
discretization parameter A and the number of kept states Nc- 
The dashed line denotes the analytic result Eq. (34) . All other 
NRG results for J in this paper were obtained at A = 4. 



where y — \/l + x"^ . The integration can be carried out 
analytically if one expands in terms of A/F, e/F or F/A, 
e/A. One obtains 



lim J = 



sin((/>) 



AJ 2v/(e/F)2+cos2(0/2)' 



(35) 



This shows that J is of order F (A) at large A (F). 
A small e 7^ affects the current mainly for cf) > 
2arccos(e/F). Running fRG, we reproduce these results. 
This is, however, only a consistency check for our numer- 
ics but not for the fRG approximation (which by con- 
struction is exact at J7 = 0) itself. In contrast, reproduc- 
ing noninteracting results is an important check for the 
NRG. Particularly for two-channel models, one cannot 
expect such a high numerical accuracy as one achieved 
for single-channel models. We have checked our NRG 
data in several ways. Fig. 6 shows an example, which is 
the Josephson current at C/ = in the particle-hole sym- 
metric case. Here, F is chosen to be much larger than 
the gap, F = 43. 2A. The NRG results are examined for 
several sets of (A, Nc) given by (□: 4.0, 700), (•: 2.0, 
700), (o: 2.0, 1716), and (x: 6.0, 700). The NRG data 
can be compared with the exact U = result (dashed 
line). As we discussed in Sec. IIIB, the numerical accu- 
racy is controlled by the discretization parameter A and 
the number of the low-energy states Nc retained in the 
truncation procedure for constructing the Hilbert space 
for the next NRG step. The results show that for A > 4, 
the NRG iterations with 700 kept states generate cor- 
rect convergent data (see the points corresponding to □ 
and X in the inset). On the other hand, the data for 
a small discretization A = 2 shows that 700 states are 
not enough, and the results (•) deviate from the exact 
ones. However, with a large number of the kept states 
[Nc = 1716), the numerical accuracy can be improved 



nicely, as the results plotted with o approach to the cor- 
rect ones. The dependence of the convergence on the 
truncation procedure becomes sensitive particularly at 
F 3> A, when the wave function of the Andreev bound 
state penetrates deep into the superconducting leads. We 
have confirmed that all statements from above also found 
for finite U. From these checks, we see that for the two- 
channel model the truncation has to be performed with 
special care, in particular for a small discretization such 
as A = 2. We have also confirmed that for larger dis- 
cretizations A ~ 4, 700 states are enough to obtain con- 
vergent results in the interacting case. Therefore, in this 
paper we have chosen A = 4 for all other NRG results for 
J, and have benchmarked the convergence at some data 
points on each of the curves by carrying out cross-checks 
using iterations with 1716 kept states. 



B. The atomic limit 

It is instructive to investigate the so-called atomic limit 
A = cx) at finite U > 0. Even though the current in 
the doublet phase vanishes at A = 00, one can obtain 
analytical expressions for J in the singlet phase as well 
as for the phase boundary. Here, we focus on symmetric 
couplings, the more general situation with F^ 7^ F^ is 
discussed in Ref. 29. 

At A = 00, the inverse of the free propagator (18) 
reduces to 



iuj ~ e Ad 
Kd ioj + e 



(36) 



where A^ = T cos{(j)/2). Including the interacting part, 
the problem becomes equivalent to solving an effective 
two-level Hamiltonian 



A 



V2 
t 1 



(37) 



Diagonalizing this operator yields a non-degenerate [dou- 
bly degenerate] ground state if B{(j>) > U/2 [B{(j>) < 



[7/2], where i3((/>) = ye^-l-A^. Hence, the boundary 
between the singlet and doublet phase is described by 



m)'-©'-cosV/2) = 0. 



(38) 



This shows that at small U the system never exhibits a 
phase transition no matter how large A. The Kondo ef- 
fect is not active and the ground state always remains a 
(BCS) singlet. As discussed above, this singlet phase sta- 
bilizes with increasing e and decreasing </>. A non-trivial 
test for the fRG approximation is to compare critical lines 
obtained numerically at large A/F with the exact result 
Eq. (38). This is done in Fig. 7(a), showing that fRG 
captures the essential behavior of the critical lines. 
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u/r u/7ir, <\>/-K 

FIG. 7: (Color online) (a) Comparison of fRG (solid) with the exact result [dashed, from Eq. (38)] for the critical phase 
difference at e = and the critical gate voltage at (fj/ir — 0.5 describing the phase boundary in the atomic limit. fRG 
results were obtained at A/F = 1000. (b) The same, but comparing the U - and <j) - dependence (at e/F = 0.5, ^/tt = 0.5 and 
e/F = 0.5, U/r — 1, respectively) of the self-energy components. Only the results in the singlet phase are shown. The end of 
the lines S(i7) and SA(t/) indicate the transition into the doublet phase. 



The Josephson current in the singlet phase follows from 
differentiating the energy of the corresponding state E = 
U /A — B{(t)) w.r.t. the phase difference (j) [see Eq. (C3)]: 



sin(0) 



2 v/(e/r)2+cos2(,/)/2) 



(39) 



This coincides with the noninteracting result Eq. (35). 
Since the energy of the doublet state E = —U /A is inde- 
pendent of 0, the current in the doublet phase vanishes. 
In Fig. 8 we show that the current (in the singlet phase) 
obtained from fRG at large A/T indeed falls onto the 
curve described by Eq. (39). 

The knowledge of the exact eigenstates also allows us 
to compute the exact Green function (using the spectral 
representation) and from this the exact self-energy. In 
the singlet phase the components arc given by 



Ue 

m<t>) 



Sa(«w) = 



(40) 



The self-energy is purely of first order in U and frequency- 
independent. In Fig. 7(b) we compare the fRG approx- 
imated self-energy to this exact result. Overall the fRG 
reproduces the parameter dependencies quite well. How- 
ever, compared to the exact result the fRG self-energy 
contains higher-order corrections which for larger U /T 
lead to deviations from the strictly linear U - dependence. 
In the doublet phase the exact self-energy is given by 



Yi(iuj) = — 



iuj + e 



4 u;2 + [S(0)]2 



Ad 



4 Lu'^ + [B{(t))]^' 
(41) 

Remarkably it is purely of second order in U and in con- 
trast to the self-energy in the singlet phase frequency- 
dependent. One cannot expect these properties to be 
reproduced by the truncated fRG which neither keeps 



all terms of order U"^ nor any frequency dependence. 
Indeed, the fRG self-energy in the doublet phase is al- 
ways (i.e. also in the atomic limit) given by S = — e and 
Sa = Ad. 

To conclude, we have demonstrated fRG to well re- 
produce analytical results for A = oo. As mentioned in 
Sec. IIIB, NRG is not as accurate for two-channel mod- 
els as it is for single-channel models, and in our case only 
the first several digits of the expectation values can be 
trusted. Thus, for calculating the very small current J/ Jc 
(of order < 10"'^) accurately in a doublet state (particu- 
larly at r ^ A for large U), it is still not efficient to use 
NRG at present. 



C. NRG vs. QMC 

The Josephson current was previously computed us- 
ing NRG by Choi et al. (CLKB, Ref. 38). The accu- 
racy of their results was questioned by Siano and Eg- 
ger (SE, Ref. 36) who performed quantum Monte Carlo 
(QMC) calculations. However, this issue has not been 
finally resolved yet.'^^''"' Furthermore, QMC is an inher- 
ently finitc-T method, whereas CLKB's NRG approach 
mainly focused on the zero temperature limit. Since both 
methods are generally regarded to benefit from numerical 
exactness, more clarifying work is needed. 

In the present paper we have re-examined the NRG 
calculation for the Josephson current checking the nu- 
merical accuracy carefully as discussed in Sec. VA. The 
parameters in Fig. 4(a) are chosen for direct comparison 
to Fig. 3 of CLKB's paper.'^6,48 whereas the data we ob- 
tained from fRG and NRG agree fairly well, there is a 
sizable discrepancy to CLKB's results. Most important, 
the amplitude of J exceeds the one computed by CLKB 
by a factor of about 1.5 - 2. Our results suggest that 
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FIG. 8: (Color online) Josephson current obtained from fRG 
(solid lines) in the atomic limit A = oo at T = 0. As fRG 
calculations were performed at A/F = 1000 < oo, the cur- 
rent in the doublet phase is finite due to F/A corrections [see 
Eq. (31)]. Note that in this phase J (which was scaled up 
by a factor of A/F) is independent of U and e. The dashed 
lines (which are partly covered by the solid lines on this scale) 
show the analytical result for the singlet phase [Eq. (39)]. 



CLKB's data for the Josephson current have a problem 
in the amplitude, although their results capture the 4> de- 
pendence correctly except in the region near (/) = tt.^^ In 
addition, we have studied the <f) - dependence of the ener- 
gies Ei of the first and second many-body excited states 
(see Fig. 9), which emerge below the gap (0 < < A). 
These in-gap excitations determine the temperature de- 
pendence at T ^ A (see below). We observe that our 
results of the first excited energy agree reasonably well 
with the peak position of the single-particle spectrum 
given in Fig. 2(d) of Ref. 38, although the broadening of 
the bound-state peak seen in the figure of CLKB must 
be artificial. One should note that the first and second 
excited states are degenerate ai (f) = tt. This is caused by 
a special symmetry holding at (/> = tt between cigenstates 
with the quantum number Q and those with —Q?^ 

In order to further clarify the discrepancies between 
the previous NRG and QMC results, we have re- 
examined the temperature dependence of the Josephson 
current using our NRG code. The results are shown in 
Fig. 10, where the parameters are taken to be the same as 
those used for the T = results at A/F = 0.37 given in 
Fig. 4(a). Thus, our data can be directly compared with 
Fig. 1 of Ref. 39. The general features of the T - depen- 
dence of CLKB's results are consistent with our findings, 
although the amplitudes again differ approximately by a 
factor of the order of 2. 

Figure 10 can also be compared with SB's QMC re- 
sults, namely with Fig. 1 of Ref. 39. We observe that 
the amplitude of the Josephson current obtained by SB 
is consistent with ours in the region 7r/2 < </> < tt where 
the ground state is a doublet. However, SB's results of 
the current are very small at < < 0.27r compared 




FIG. 9: (Color online) Energies Ei of the first and second 
many-body excited state emerging below the edge of the su- 
perconducting gap at U /V — 5.2, and A/F = 0.37. Ei is 
measured with respect to the ground-state energy. {Q, S) de- 
notes the set of quantum numbers introduced in Sec. IIIB. 
The NRG parameters are taken to be same as in Fig. 4(a). 



to our NRG data. The reason for the discrepancy may 
be inferred from the (j> - dependence of the excitation en- 
ergy. As shown in Fig. 9, the first excitation energy Ei 
is very small (0 < £^1 < O.IA) for < < 7r/2. In 
particular, Ei is smaller than (or almost equal to) the 
temperature T = O.IA that SB used for their calcula- 
tions. They have set up the transition probability for the 
importance sampling introducing the ultraviolet cutoff 
for the summations over the Matsubara frequencies such 
that \uJn\ < ttPT, where P is the Trotter number. ^^.so 
Therefore it is questionable whether the transition prob- 
ability captures accurately the information about the in- 
gap states if the excitation energies are very small. This 
will not cause any problems if Ei is larger than T, and it 
explains naturally the agreement between the QMC and 
our NRG data for 7r/2 < (/) < tt. 



D. Mean- field 

For the ordinary single impurity Anderson model it is 
well-known that the Kondo effect cannot be described 
within a mean-field framework. Despite this fact the 
Hartree-Fock approximation was used to compute the 
Josephson current through a single impurity coupled to 
BCS leads. However, previous approaches were either in- 
complete (the induced anomalous term was discarded in 
Ref. 31) or inaccurate (Ref. 32; see below). For reference, 
we have thus performed our own Hartree-Fock calcula- 
tion based on self-consistent equations which we derived 
in agreement with Yoshioka and Ohashi (YO, Ref. 32). 
As within truncated fRG, the Hartree-Fock self-energy 
is frequency- independent. Astonishingly, the general pic- 
ture of the phase transition is captured on this mean-field 
level, but the quantitative agreement with reliable NRG 
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FIG. 10: (Color online) NRG results for the Josephson current 
at finite temperatures T. The parameters are the same as in 
Fig. 9. 



FIG. 11: (Color online) The same as in Fig. 4(b), but obtained 
within a self-consistent Hartree-Fock framework. 



data is poor [compare Figs. 4(b) and 11]. However, the 
observation that Hartree-Fock succeeds in quaUtativcly 
describing the behaviour of the Josephson current is an 
uncontrolled result since an approximation which inher- 
ently does not contain the Kondo temperature cannot 
be expected to describe a transition governed by A/Tk- 
In addition, the appearance of a phase with J < is 
caused by an unphysical breaking of the spin symme- 
try (the ground state is not a doublet). Thus, a direct 
comparison between Hartree-Fock and NRG is actually 
of limited value and provided by Fig. 11 only for reasons 
of completeness. 

It is important to mention that even though we derive 
the same mean- field equations as YO, we cannot repro- 
duce their numerical solution. In particular, YO observe 
the induced magnetization (the difference between the 
average occupation numbers of spin up and down elec- 
trons) to increase continuously from zero when U is in- 
creased. In contrast, our solution shows a discontinuous 
jump when the system enters the phase with J < 0. We 
have double-checked our data using different routines to 
solve the self-consistency problem and arc thus tempted 
to conclude that YO's results arc inaccurate due to nu- 
merical issues. 



VI. CONCLUSIONS 

In this paper we have presented phase diagrams and 
the Josephson current for a single Anderson impurity 
coupled to BCS superconducting leads. We obtained our 
data using the frameworks of the functional and numeri- 
cal renormalization group, respectively. The well-known 
accuracy of the NRG was established by comparisons 
with the exact solution at ?7 = 0, allowing us to show that 
previous NRG results for the Josephson current by Choi 
et al. are smaller approximately by a factor 2 compared to 
the correct values. Using NRG wc have also re-examined 



the temperature dependence of the current. Our results 
agree reasonably well with the QMC data by Siano and 
Egger when the in-gap excited energy E is larger than 
T, while in the opposite case < < T < A the QMC 
results become less accurate. We have demonstrated the 
truncated fRG to perform well compared with NRG at 
small to intermediate interactions and to reproduce an- 
alytical results which we derived in the limit A = oo. 
As fRG requires virtually no numerical resources it can 
be used to fast provide current - phase curves needed to 
extract physical quantities out of experimental data.^^ 
Concerning this issue it would be desirable to extend the 
fRG scheme to treat finite temperatures. Doing this in a 
reliable but at the same time numerically efficient way is 
difficult and needs further investigation. 



Acknowledgments 

This work was initiated by enlightening discussions 
with T. Novotny and J. Paaske when V.M. was visit- 
ing the Nano-Science Center of the University of Copen- 
hagen. V.M. very much enjoyed the hospitality of this 
institution. A.O. is grateful to the condensed matter 
physics group of the Universitat Gottingen for the kind 
hospitality during his stay. His visit was supported by 
the Deutsche Forschungsgemcinschaft (SFB 602). We are 
much obliged to S. Andcrgassen, J. Bauer, A. C. Hcwson, 
Y. Nisikawa, T. Novotny, T. Pruschke, K. Schonhammer, 
Yoichi Tanaka, and Yoshihide Tanaka for valuable dis- 
cussions and thank P. Plotz for careful reading of the 
manuscript. C.K. and V.M. are grateful to the Deutsche 
Forschungsgemcinschaft (FOR 723) for support. A.O. is 
supported by the Grant-in-Aid for Scientific Research 
from JSPS. NRG calculations were partly carried out on 
SX8 at the computation center of the Yukawa Institute. 



13 



APPENDIX A: NRG APPROACH 

In the NRG approach, the continuous conduction 
bands of the original model Hamiltonian are discretized 
and described by the fermions „^ and defined 
on a linear chain for n > on the left and right, 
respectively.^^ In the case of equal gaps (A^ = Aj^ = A) 
and symmetric couplings (F^ = F/j = r/2) it is conve- 
nient to use the linear combinations'^* 



i0/4 J" 



^/2 

J H^ncr ^ J L,na 



(Al) 



-^/2^ 

The discretized Hamiltonian of the NRG, 



TT _ rrdot , frlcad , rrcoup 

^NRG - tL + i^NRG + ■'^NRG' 



(A2) 



can be expressed in a way such that it contains only real 
parameters for arbitrary phase difference (j), even though 
the operator H defined in Sec. II is generally a complex 
Hermitian matrix when <f> is finite. Namely, the part of 
the Hamiltonian modeling the BCS leads can be written 



lead 
NRG 



E 



n— L cr 



(A3) 



The hopping matrix elements fall off exponentially with 
the distance n from the dot: 



Un = D 



1 + 1/ A (1-1/A"+^)A 



(A4) 



2 ^1 - 1/A2»+Vl - l/A2«+3 
The discretized version of the coupling Hamiltonian reads 

-f^NRG = V2UNRG COs(0/4) ^ (sj^d^ + rfj, SOo 

(T 

+%/2uNRG sin(0/4) ^ (plada + dlpoa j , 

(A5) 



where we have defined wnrg = v^2_DFnrg/^i and 
^aF 



NRG 



11 + 1/A, ^ 

Aa = '— In A. 

21-1/A 



(A6) 



The factor A\ is necessary for correctly reproducing the 
original model in the continuum limit A 24,25,52 rpj^g 
discretized version of the current operator corresponds 
to a derivative of the timneling Hamiltonian w.r.t. the 
phase difference 0: 



Jnrg = 290i/^'^Q. 



(A7) 



For a numerical implementation it is useful to exploit 
symmetry properties in order to reduce the size of the 
matrices to be diagonalized in each NRG step. Particu- 
larly, the pseudo-spin Q^^ defined by^^ 



d\dl 



n=0 



H.C. 



(A8) 

is conserved, [QP'',i?NRG] = 0, in the particle-hole sym- 
metric case (e = 0). Therefore, the eigenvalue Q of the 
operator Qp^ can be used as a quantum number to clas- 
sify the Hilbert space in addition to the total spin S asso- 
ciated with the rotational symmetry of the real spin.^'^'^^ 



APPENDIX B: CURRENT FORMULA 

In this section we derive an exact formula relating the 
Josephson current to the self-energy. To this end, we 
define a current operator at lead s = L, R as usual. 



Js^dtN,=i[H,Ns 



(Bl) 



Two terms of the Hamiltonian fail to commute with the 
particle number operator Ns'. 



H.C. 



(B2) 



The expectation value of the second term vanishes since 
Ase"^= ^ J^ki^^kiCs-ki)- The first term can be evalu- 
ated using a projection technique for the noninteracting 
Green function at the dot-lead interface:^'^ 



(B3) 



with = 1, 2 denoting Nambu indices of the dot and 

of the local site at the end of lead s, respectively. The 
generalization of Eq. (B3) to the interacting Green func- 
tion is achieved straight-forwardly by virtue of the Dyson 
equation: 



(B4) 



Using Eqs. (B3) and (B4) to evaluate the expectation 
value of Eq. (B2) we obtain the following expression for 
the current (see also Ref. 29): 



{Js 



At 



r [^2,i(i^)gb(* 



(B5) 



with g{z) being the dot Green function, and (3 = 1/T. 
The most general form (in absence of a magnetic field) 
of g{iu)) is 



g{iuj) 



iu) + e + Y,{iuj)* Y,j\{iLo) — IS. 



D{itjj) \ Ti/\{iLo)* — A* iuj ^ e — Yj{iuj) J 

(B6) 



14 



where D{iuj) denotes the determinant 



(B7) 



The self-energy components fulfiU the symmetry relations 
= Yi{iuj)* and Sa(— iw) = SA(*ti^),'^'' implying 
that D{iuj) is purely real. Employing Eqs. (B5) and (B6) 
allows for recasting the current formula into the simple 
form 



{Js 



r,A, Im [c-^-^^EaM] 1 



(B8) 



where we have introduced the notation L = R,R = L. 
One should note that this is an exact result, the gener- 
alization to the case of broken spin symmetry (i.e. in 
presence of a magnetic field) being achieved straight- 
forwardly. At zero temperature, the Matsubara sum 
can be evaluated as an integral by replacing 1//3 ^ — > 
l/27r J duj. Within fRG and for symmetric parameters, 
Sa is real and we obtain Eq. (27). 



APPENDIX C: CURRENT CONSERVATION 

This section is devoted to the question of current con- 
servation, an issue which can also be tackled within a 
more general framework using generating functionals.^^ 
Within the model Hamiltonian Eq. (1), electrons cannot 
be created or annihilated on the quantum dot. Hence, 
one would expect the Josephson current to be conserved: 



{Jl) = -{Jr). 



(CI) 



This can indeed be shown analytically by applying 



a gauge transformation Cgka 
(.-i't'R/-ifla, and H ^ H, with 



^sko 



(C2) 



The current operator Eq. (B2) can then be expressed as 
a derivative of the grand canonical potential f2 w.r.t. the 
phase difference = 0l — 0^: 



{Jl) = 2{d^H{cf,)) = 2d^n{cl>). 



(C3) 



The same result is obtained for — (J_r). Thus, the Joseph- 
son current is conserved. 

Current conservation implies a symmetry relation for 
the self-energy. In particular, plugging Eq. (B8) into 
Eq. (CI) yields 



Im 



E E 

s—L^R iuj 



D{iuj) 



0. 



(C4) 



This equation is fulfilled if the interacting many-particle 
system is solved exactly. On the other hand, any ap- 
proximate method to calculate the self-energy is current- 
conserving if and only if Eq. (C4) holds. 

To discuss the issue of current conservation for 
frequency-independent approximations (such as trun- 
cated fRG and Hartree-Fock), it is instructive to rewrite 
the off-diagonal component of the self-energy in terms of 
a function g{iuj) ~ g{—iLu) defined by 



^Aiioj) = g{ioj) ^ ^ as{i 



(C5) 



where as{iuj) = l^a{i^) / D{iuj). Equation (C4) can then 
be recast as 

^™ E E «^i(«^)«s2(*f')5*(*^) = 0- (C6) 

It follows that for any frequency-independent approxima- 
tion to the self-energy, current conservation is equivalent 
to 



Im g = 0. 



(C7) 



It is easy to show that this condition is always fulfilled if 
the problem is treated within a self-consistent Hartree- 
Fock approach. In contrast, the Josephson current is 
only conserved for symmetric gaps (A^, = An) if the self- 
energy is computed from the fRG flow equations (23) and 
(24). This can be seen by using g^ = / as{ii') 
to describe the flow of the off-diagonal part of the self- 
energy. The flow equation is determined by Eq. (24), 



dAg' 



9 



(C8) 



The second term on the r.h.s. is real only for A^ = Aij. 
Hence, within the fRG framework the current is con- 
served if and only if the superconducting gaps of the left 
and right lead are equal. 

Further insight into the structure of the self-energy can 
be gained by requiring Eq. (C4) to hold term by term. 
This is only possible if 

Sa(«w) = fiiuj)A{iuj), f{iuj) = fi-iuj) e M. (C9) 

Thus, it is reasonable to consider the flow of = 
S^/^^Fse"^" instead of E^. For symmetric gaps 
(Al = Ar = A), Eq. (24) gives 



A 



VA2 + A2 



(CIO) 



Since the r.h.s. is real so is the function f^, and E^ is of 
the form (C9). Hence, for symmetric SC gaps the identity 
(C4) is fulfilled term by term within the fRG approach. 
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